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Abstract 

A  current  thrust  in  medical  research  is  the  development  of  a  non-invasive  method  for  detection,  localization, 
and  characterization  of  an  arterial  stenosis  (a  blockage  or  partial  blockage  in  an  artery).  A  method  has  been 
proposed  to  detect  shear  waves  in  the  chest  cavity  which  have  been  generated  by  disturbances  in  the  blood  flow 
resulting  from  a  stenosis.  In  order  to  develop  this  methodology  further,  we  use  one-dimensional  shear  wave 
experimental  data  from  novel  acoustic  phantoms  to  validate  a  corresponding  viscoelastic  mathematical  model. 
We  estimate  model  parameters  which  give  a  good  fit  (in  a  sense  to  be  precisely  defined)  to  the  experimental  data, 
and  use  asymptotic  error  theory  to  provide  confidence  intervals  for  parameter  estimates.  Finally,  since  a  robust 
error  model  is  necessary  for  accurate  parameter  estimates  and  confidence  analysis,  we  include  a  comparison  of 
absolute  and  relative  models  for  measurement  error. 


Mathematics  Subject  Classification:  62F12;  62F40;  65M32;  74D05. 

Key  words:  viscoelastic  model;  sensitivity  analysis;  inverse  problem;  asymptotic  theory. 
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1  Introduction 


Coronary  artery  disease  (CAD)  is  an  increasingly  prevalent  medical  condition,  often  a  precursor  to  and  cause  of 
a  patient  experiencing  cardiac  arrest.  Current  methods  for  detection  of  arterial  stenoses  (blocked  arteries)  include 
the  angiogram  and  CT  scans.  Angiograms  are  viable  but  quite  invasive,  while  CT  scans  are  expensive,  introduce 
radiation  into  the  patient,  and  can  only  detect  hard  plaques  (blockages).  A  desirable  new  detection  method  would 
be  noninvasive  and  less  expensive  but  still  effective.  To  this  end,  using  acoustic  waves  generated  by  stenoses  has 
been  proposed.  This  would  place  sensors  on  the  surface  of  the  chest  to  listen  for  sounds  from  coronary  arteries,  with 
the  hope  of  detecting  and  then  localizing  any  blockages. 

The  current  understanding  (see,  e.g.,  [8,  26,  45])  of  the  process  is  that  turbulent  flow  produces  normal  forces 
on  the  vessel  walls  at  and  downstream  of  a  stenosis,  which  then  exert  pressure  on  the  vessels  wall  causing  a  small 
displacement  in  the  surrounding  soft  tissue.  Previous  work  (e.g.  [2,  32,  46,  43,  44,  52,  56])  has  demonstrated  the 
existence  of  such  sounds,  and  that  they  are  possible  to  detect  on  the  surface  of  the  chest.  We  thus  see  that  the 
system  couples  two  processes:  (1)  the  generation  of  pressure  and  shear  waves  transmitted  into  the  body  through  the 
arterial  wall  as  a  result  of  the  turbulent  blood  flow  generated  by  a  stenosis,  and  (2)  the  propagation  of  pressure  and 
shear  waves  through  the  chest  to  sensors  attached  to  the  chest  wall.  The  first  process  is  not  completely  understood, 
though  some  ideas  are  present  in  the  literature.  Various  researchers  [24,  25,  26,  36,  39,  41,  45,  48,  49,  55,  59]  have 
directly  examined  modeling  blood  flow  through  arteries.  All  have  attempted  to  characterize  the  turbulence  in  the 
flow,  which  some  then  used  to  examine  the  sound  field  propagated  into  the  chest.  In  the  current  work,  we  will 
not  focus  on  aspects  of  turbulent  flow,  leaving  that  instead  as  an  input  to  be  later  properly  characterized  when  the 
stenosed  artery  situation  is  modeled  directly.  We  will  focus  on  the  second  process,  understanding  the  propagation  of 
sounds  through  the  chest  cavity  which  result  from  stenosed  coronary  arteries. 

The  modeling  and  detection  of  waves  transmitted  through  the  body  has  been  approached  in  different  ways.  One 
approach  has  focused  on  characterizing  properties  of  the  sounds  detected  on  the  surface  of  the  chest,  characterizing 
aspects  such  as  primary  frequencies  that  these  sounds  exhibit.  This  line  of  work  has  been  studied  by  Semmlow,  et 
ah,  [2,  3,  4,  5,  6,  47,  53,  54,  61],  as  well  as  by  other  groups  [24,  25,  26,  34,  46,  56,  57,  58].  Their  methods  are  based 
on  general  sound  features  and  detection  through  statistical  methods,  rather  than  modeling  the  underlying  physics  of 
sound  transmission  through  the  body.  These  methods  have  the  benefit  of  being  fast  and  fairly  simple  to  implement, 
but  do  not  provide  a  characterization  for  the  mechanisms  of  wave  transmission. 

In  another  direction,  more  relevant  to  the  situation  this  paper  will  study,  researchers  have  worked  to  model 
the  physics  of  sound  wave  transmission  through  the  body.  As  is  common  in  many  physical  wave  phenomena,  both 
pressure  and  shear  waves  propagate  into  the  body  as  a  result  of  stenosed  coronary  arteries.  Since  shear  waves  in 
general  have  a  lower  amplitude  than  pressure  waves,  intuition  might  suggest  detecting  pressure  waves  should  be  our 
focus;  past  research  indicates  the  opposite  is  true  when  seeking  to  detect  stenoses.  Various  groups  [28,  30,  40,  45] 
have  demonstrated  that  shear  waves  should  be  the  focus  of  detection  efforts.  The  frequency  ranges  for  shear  waves 
resulting  from  coronary  stenoses  are  below  2000  Hz  [25,  28,  29,  33,  51,  53].  In  the  range  of  these  frequencies,  the 
pressure  waves  propagate  very  quickly  while  shear  waves  propagate  much  more  slowly,  which  in  practice  means  that 
pressure  waves  are  difficult  to  measure.  Furthermore,  the  lower  speed  of  shear  waves  implies  a  shorter  wavelength 
for  a  given  frequency,  thus  providing  better  spatial  resolution  when  attempting  to  localize  a  stenosis.  Also,  in  the 
context  of  waves  propagating  in  tissue  or  tissue- mimicking  materials,  shear  waves  are  measurable  at  greater  distances 
from  the  source  of  the  disturbance  (see,  e.g.,  Figures  9b  and  10  of  [28]).  Thus,  in  this  work  we  also  focus  on  studying 
a  model  for  shear  wave  propagation  (though  results  where  we  also  test  a  pressure  wave  model  are  available  in  a 
technical  report  [14]). 

The  benefits  of  using  a  viscoelastic  wave  propagation  model  in  various  contexts  have  been  previously  studied 
[28,  30,  31,  33,  37,  40,  51].  In  these  references,  the  authors  focused  on  determining  the  elastic  modulus  and  viscoelastic 
parameter  based  on  the  shear  wave  speed  and  attenuation  in  either  a  gel  mold  or  physical  tissue,  in  both  a  stenosis 
context  and  general  tissue  shear  wave  propagation.  The  models  were  developed  using  plane  waves  in  such  a  way 
that  algebraic  expressions  were  available  for  shear  wave  speed  as  a  function  of  frequency,  elastic  modulus,  and  the 
viscoelastic  parameter.  These  demonstrate  that  modeling  the  underlying  physics  is  not  only  possible,  but  quite 
beneficial  in  understanding  shear  wave  propagation.  These  investigators  also  showed  that  a  Kelvin- Voigt  damping 
model  is  most  appropriate  for  the  situation;  we  will  incorporate  this  into  our  model  (and  will  discuss  it  more  fully 
below). 

In  this  work,  we  will  take  the  physical  models  further,  developing  a  dynamic  model  of  the  shear  waves  propagating 
through  a  tissue-mimicking  material.  Our  goal  will  be  to  use  this  viscoelastic  model  with  data  from  a  tissue- mimicking 
homogeneous  gel  mold  to  validate  the  model  and  understand  the  uncertainty  in  the  model  parameters.  The  model 
here  will  incorporate  the  standard  elastic  modulus  and  bulk  viscoelastic  parameter,  as  well  as  internal  variables 
governed  by  relaxation  times  which  can  be  used  to  model  how  different  portions  of  the  medium  relax  in  different 
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ways  from  being  stressed.  Wc  will  develop  the  model  in  such  a  way  as  to  allow  for  multiple  internal  variables.  In 
future  studies  these  can  be  used  to  model  the  bulk  effects  of  different  types  of  tissue,  which  would  be  closer  to  the 
in  vivo  detection  problem. 

The  model  here  will  continue  from  a  previous  line  of  work  by  Banks,  et  al.,  [8,  18,  19,  20,  21,  42,  50].  These 
models  allow  for  a  characterization  of  shear  waves  resulting  from  coronary  stenoses,  which  will  assist  in  uncovering 
the  coronary  artery  sounds  from  the  noisy  background  in  the  body.  Initial  experiments  were  conducted  where  a 
gel  mold  was  built  with  a  tube  running  through  the  middle;  cases  where  the  tube  was  unblocked  were  compared  to 
those  with  partial  blockages,  and  the  results  suggested  that  there  were  significant  differences  in  sound  generation 
between  the  blocked  and  unblocked  cases.  Unfortunately,  this  line  of  work  ended  before  experimental  data  could  be 
incorporated  and  used  to  validate  models.  The  current  work  returns  to  a  one-dimensional  model  and  experimental 
setup,  developing  a  model  closely  related  to  that  in  [8]  which  we  validate  with  lab  data  from  a  homogeneous  gel 
phantom  and  use  to  examine  uncertainty  in  our  parameter  estimates.  This  work  will  examine  confidence  in  the 
relaxation  times,  which  will  demonstrate  that  the  addition  of  internal  variables  into  the  model  is  viable. 

We  now  summarize  the  goals  and  direction  for  the  current  work.  Here,  we  incorporate  all  of  the  aforementioned 
physical  modeling  concerns  and  continue  the  work  of  our  concept  paper  [13]  by  focusing  on  wave  propagation  through 
a  homogeneous  viscoelastic  medium.  The  constitutive  relationship  will  be  slightly  more  general  than  in  our  concept 
paper,  and  is  based  on  Fung’s  viscoelastic  formulation  that  has  been  validated  by  actual  tissue  experimental  data 
(more  details  below).  This  relationship  will  be  used  in  one-dimensional  shear  wave  dynamical  models  to  compute 
inverse  problem  results  for  the  one-dimensional  case  using  experimental  data  from  a  tissue-mimicking  gel  mold.  This 
data  comes  from  novel  acoustic  phantoms  built  and  tested  at  Queen  Mary,  University  of  London  (QMUL)  and  Barts 
Health  Trust  (BHT)  in  England.  The  data  will  be  taken  in  such  a  way  that  measured  displacements  upon  releasing 
a  load  are  primarily  affected  by  the  properties  of  the  medium.  Displacement  data  for  our  case  encodes  the  frequency 
and  attenuation  properties  of  the  medium,  which  we  then  quantify  through  model  parameter  estimation.  This  will 
firmly  set  the  groundwork  for  future  investigations  into  validating  the  model  in  a  more  complex  medium,  as  well 
as  investigation  into  a  proper  description  of  an  input  for  the  wave  propagation  model;  this  input  would  result  from 
models  of  the  sounds  (arterial  wall  displacements)  due  to  turbulent  blood  flow  resulting  from  a  stenosis. 

2  Experimental  Setup 

A  novel  experiment  has  been  devised  at  QMUL  to  gather  one-dimensional  shear  data.  Devices  have  been  designed 
(see  left  pane  of  Figure  1),  in  which  an  agar  gel  mold  phantom  (homogeneous,  97%  water,  density  p  =  1010  kg/m3)  is 
loaded  into  the  rig,  a  weight  is  attached  applying  stress  to  the  phantom,  and  then  the  weight  is  released,  causing  the 
material  to  oscillate.  The  displacement  motion  of  the  material  throughout  the  experiment  is  measured  with  a  laser 
device.  This  motion  is  due  to  the  material  response  to  the  weight  release,  and  thus  encodes  information  about  the 
material  characteristics  (elastic  modulus,  damping  characteristics,  etc.).  The  choice  of  loading  and  a  quick  release  is 
designed  to  produce  dynamic  data;  the  idea  was  inspired  in  part  by  the  impacts  the  stenosed  vessel  wall  experiences 
with  each  heartbeat  and  also  by  past  success  in  gathering  shear  data  for  filled  rubber  elastomers  using  an  initially 
loaded  rubber  sample  which  then  underwent  an  impulsive  hammer  hit  (see  e.g.  [18,  19]).  This  yields  one-dimensional 
shear  data  in  the  radial  direction  perpendicular  to  the  vertical  axis  (right  pane  of  Figure  1). 


Figure  1:  Shear  configuration,  where  LDT  denotes  the  laser  displacement  transducer,  (left)  Experimental  setup  of 
agar  phantom,  (right)  Schematic  with  one-dimensional  domain  denoted. 
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In  this  work,  we  focus  on  results  generated  from  a  load  weight  of  264  g,  as  this  weight  level  produces  a  well-defined 
response.  In  the  future,  we  may  also  incorporate  data  from  smaller  weight  levels  in  order  to  examine  the  different 
phantom  responses  to  other  weight  levels.  The  gel  phantoms  were  stored  in  water  when  not  in  use,  which  keeps  the 
gel  at  the  desired  97%  water  composition. 

When  the  experiment  is  conducted,  data  like  those  depicted  in  Figure  2  are  produced.  The  material  is  at  rest, 


Figure  2:  Sample  one-dimensional  data.  Loading  of  the  material  (initially  at  rest)  begins  at  t  =  Ti,  and  the  material 
is  loaded  and  continues  to  relax  for  t  £  (r2,r3).  At  time  T3  the  load  is  cut  which  takes  roughly  10-15  ms.  The  gel 
is  then  freely  oscillating  at  T4  =  0,  and  oscillations  continue  for  a  period  of  time  dependent  on  the  loading  weight. 
The  value  A  is  the  displacement  of  the  material  at  the  beginning  of  free  oscillations.  The  overall  displacement  scale 
of  the  data  is  on  the  order  of  10-4  m,  while  the  oscillations  immediately  after  the  weight  release  are  on  the  order  of 
10-5  m. 

a  weight  is  added  and  allowed  to  settle,  then  the  string  holding  the  weight  is  rapidly  cut  with  a  flame  to  allow  the 
material  to  freely  oscillate.  Once  oscillations  have  died  out  the  material  relaxes  back  toward  a  stable  state.  The  key 
pieces  that  will  be  modeled  are  the  loading  profile  (loading  begins  at  t  =  Ti  and  lasts  until  the  weight  begins  to  be 
released  at  t  =  Ts),  which  we  will  model  as  instantaneous  loading  to  position  A,  and  the  oscillations  after  weight 
release  (free  oscillations  begin  at  t  =  T4  =  0)  which  are  the  main  object  of  investigations  here.  For  more  information 
on  the  experimental  setup,  interested  readers  may  refer  to  [27]. 


3  Model  Development  and  Constitutive  Equation 


With  the  setup  of  the  experiment  in  mind,  we  can  turn  to  our  mathematical  model  of  wave  propagation.  The  model 
will  be  developed  to  take  into  account  all  features  of  the  data,  including  the  loading  profile  and  the  relaxation  present 
in  data.  Since  our  phantom  is  cylindrical,  the  model  development  begins  with  three-dimensional  equations  of  motion 
in  cylindrical  coordinates.  These  are  given  in  [50,  p.20] ,  and  also  in  [42],  and  are  derived  from  momentum  and  mass 
balance  principles.  Using  the  fact  that  the  gel  is  homogeneous  and  that  there  are  symmetries  in  the  experimental 
design,  these  three-dimensional  equations  can  be  reduced  to  simplified  one-dimensional  model. 

Let  u(r,t)  represent  the  displacement  of  the  material  at  position  r  and  time  t ,  with  r  £  (rmin,  rmax)  (for  our 
device,  rmin  =  0.0105  m  and  rmax  =  0.054  m)  and  t  >  Ti,  where  the  time  Ti  is  chosen  as  the  beginning  of  any 
stress-strain  history  in  the  material,  and  we  are  assuming  the  material  has  been  at  rest  long  enough  that  it  is  only 
affected  by  displacements  for  t  >  Ti.  Then  the  governing  partial  differential  equation  (PDE)  becomes 


d2  d  ,  a  (r,t) 

'U’ij'max  5  —  0 

u(r,Ti)  =  0,  ut(r,Ti)=0, 


(1) 


where  p  is  the  density  of  the  material,  a  denotes  the  stress,  g(t )  is  a  function  that  describes  the  loading  process 
(to  be  discussed  later).  In  order  to  complete  these  models,  we  must  provide  a  form  for  cr.  This  is  the  constitutive 
relationship,  also  called  the  stress-strain  law  since  it  relates  strain  (^)  and/or  the  strain  rate  to  stress  a.  The  next 
sections  discuss  this  aspect  of  the  model. 
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3.1  Constitutive  equation 

We  incorporate  the  previous  modeling  ideas  together  into  a  new  constitutive  equation  for  the  shear  wave  PDEs  (1). 
We  develop  the  constitutive  equation  which  will  be  used  in  this  work  assuming  that  we  will  be  solving  the  model 
starting  at  t  =  Ti  and  thus  incorporating  both  the  loading  process  and  oscillations  into  our  dynamic  equations,  for 
the  time  being.  In  Section  3.1.5,  we  will  make  an  approximation  to  the  loading  process  which  will  allow  us  to  focus 
on  the  dynamic  oscillations  of  the  material  after  the  weight  is  released,  which  is  our  true  interest. 


3.1.1  Fung’s  quasi-linear  viscoelastic  formulation 

Some  of  the  initial  investigation  into  the  viscoelastic  nature  of  tissue  was  completed  by  Fung  [35].  His  work  is  of 
particular  interest  because  it  was  validated  in  actual  tissue.  Fung  developed  a  “quasi-linear”  model 


a(t) 


as 


(2) 


with  a  kernel  of  the  form 


G{t) 


1  +  cj^2  jrexp(— t/r)dT 
1  +  c1ii(t2/ti) 


(3) 


Within  (2),  A  represents  the  stretch  of  a  material  (A  =  l  +  ur)  and  ae  describes  the  elastic  response  to  the  elongation 
A,  given  by  (see  [35,  Sec.  7.6]) 

<re(A)  =  -p+peaUr 


where  a  and  /3  are  constants  to  be  estimated.  The  parameters  t\  and  t2  are  lower  and  upper  bounds,  respectively, 
on  relaxation  times ,  which  describe  the  ways  in  which  the  material  responds  to  imposed  stresses  and  strains.  This 
model  incorporates  a  continuum  r  £  [ti,t2]  of  relaxation  times,  which  Fung  found  to  be  necessary  in  order  for  his 
model  to  match  the  response  of  tissue,  as  well  as  a  constant  term  in  the  kernel.  This  Fung  kernel  will  serve  as  a 
baseline  for  reference  when  developing  the  model  for  this  paper. 


3.1.2  Linearized  constitutive  equation 

One  could  keep  nonlinearities  in  the  constitutive  equation  (2).  However,  we  found  (as  we  shall  see  later)  that  a  linear 
constitutive  relationship  gives  a  reasonable  approximation  to  the  data  provided  by  QMUL  and  BHT.  To  that  end, 
we  will  use  the  first  two  terms  of  the  Taylor  expansion  of  eaUr  to  approximate 

<re  ss  — /?  +  /?( 1  +  aur)  =  /3aur  =  £ ur  (4) 

where  we  have  combined  £  =  f3a  into  a  single  parameter  to  be  estimated;  £  will  be  incorporated  into  other  parameters 
later  in  model  development.  We  can  then  linearize  (2)  by  using  (4),  add  a  Kelvin- Voigt  damping  term  (a  common 
linear  viscoelastic  damping  model  [12]),  and  obtain 

<j(t)  =  Giurt(t)  +  qJ  Q(t-  s^ds  (5) 

where  Q(t)  is  a  kernel  to  be  specified.  To  an  extent,  the  Kelvin- Voigt  term  describes  the  overall  nature  of  the  damping 
present  in  the  material,  while  the  kernel  Q(t)  will  incorporate  different  material  responses  at  both  the  macroscopic 
and  microscopic  levels. 


3.1.3  Existence  and  uniqueness  for  shear  wave  model 

Before  moving  on  to  the  specific  form  of  the  constitutive  equation  kernel,  we  first  establish  existence  and  uniqueness 
for  the  shear  wave  equations  (1)  with  the  constitutive  equation  (5).  To  that  end,  we  set  up  a  similar  framework  as  in 
the  concept  paper  [13]  and  connect  those  results  to  the  current  work.  We  will  require  that  the  following  assumptions 
hold: 

(Al)  The  boundary  condition  function  satisfies  g  £  £2(Ti,T); 

(A2)  The  kernel  Q  is  differentiable  with  respect  to  t  and  with  constants  c\  and  c2  such  that  |£/(f)|  <  c\  and  \G(t)\  <  c2 
for  all  t  £  [Ti,T], 
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Let  H  =  C2(rmin,  rmax),  Y  =  {<; b\(j>  G  'H1(rmir!,,  rmax),  4>(.rmax)  =  0},  and  V*  denote  the  topological  dual  space 
of  V.  We  identify  H  with  its  topological  dual  HI*  and  thus  again  obtain  V  HI  =  HP  Y*  as  a  Gelfand  triple. 
Let  Cu,(ri,T;V)  denote  the  set  of  weakly  continuous  functions  in  V  on  [ri,T],  and  Ct  =  {v  :  [Ti,  T]  — »  HI  |  v  G 
Cu,(ri,T;V)n/:2(ri,T;Y)  and  vt  G  CW(T\,  T;  H)  (~l  £2(Ti,  T;  V)}.  The  notion  of  weakly  continuous  (i.e.,  um  — >  u  in 
Cw(Pi,T;  V))  means  that  um  — >  u  weakly  in  V  and  uniformly  in  t  G  [ri,T].  Then  a  weak  solution  u  G  Ct  for  the 
shear  equation  must  satisfy 


0  =  p(ut(t),r)t(t))  -  p  (uB(s),r]s(s))ds  +  g(s)r)(rmin,  s)ds  +  Gi  /  (usr{s),  rjr(s))ds 

Jr1  J  ri  .  J  ri 

frma x  urt{r,s) 


+C  Jr  \Jr  ^s-^Ur{mMs))dS~Gl  Jr  J 

_</r.  II  (/r.  I  “ 


-r](r,  s)drds 


(6) 


for  any  t  G  [Ti,  T]  and  rj  G  £t  and  where  (-,  •)  is  the  usual  inner  product.  Since  rmin  >  0,  there  are  no  singularities  in 
the  final  term  in  (6),  and  the  kernel  integral  in  the  numerator  of  that  term  will  converge  in  the  same  manner  as  the 
preceding  kernel  integral.  Thus,  the  arguments  from  [13]  apply  in  our  case  here,  and  we  have  the  following  theorem: 


Theorem  3.1.  Assuming  (Al)  and  (A2),  the  shear  equation  (1)  with  constitutive  relationship  (5)  has  a  unique  weak 
solution  on  any  finite  interval  [Ti,T]. 


3.1.4  Form  for  constitutive  equation  kernel  Q(t) 

We  will  now  state  the  particular  kernel  used  for  this  current  work,  and  then  manipulate  it  into  a  form  that  gives 
more  physical  insight  and  which  will  later  allow  for  a  conceptual  framework  using  internal  variables.  We  develop 
this  kernel  from  a  different  perspective  than  that  given  in  [13],  but  the  resulting  form  will  be  quite  similar.  Using 
the  notation  and  parameter  conventions  of  [12],  we  define  the  kernel  in  this  work  to  be 

Q(t-P)  =  Kr  +  K{t]P)  (7) 

where  Kr  is  a  positive  constant  representing  an  instantaneous  relaxation  modulus  (justified  by  the  fact  that  our  gel 
phantom  acts  partly  as  a  solid)  and  K(t;P)  =  J’7-exp(— t/r)dP(r)  represents  a  continuum  of  distributed  relaxation 
times  with  T  =  [ti,T2]  C  (0,  oo)  and  where  P(t)  is  a  probability  measure  on  T.  Note  that  this  form  for  Q  satisfies 
\G(t)\  <  ci  with  Q  clearly  differentiable  and  \G(t)\  <  C2  for  some  constants  ci,  C2  so  that  assumption  (A2)  is  satisfied. 
It  is  also  worth  noting  here  that  our  proposed  kernel  form  (7)  is  similar  to  that  in  Fung’s  model  (3),  as  we  see  that  Kr 

serves  as  an  analog  to  the  constant  portion  of  Fung’s  kernel  (i.e.,  ^  ^  ^ — -j — y)  and  the  K(t\  P)  portion  is  similar 

to  the  the  continuous  relaxation  spectrum  in  Fung’s  model  (i.e.,  ^ Tl  T — - — ^ ] — ). 

1  +  cln(T2/Ti) 

We  substitute  (7)  into  (5)  and  manipulate  the  form  of  the  stress  with  integration  by  parts,  noting  that  wr(ri)  =  0 
since  the  material  is  initially  at  rest  and  using  the  fact  that  AT(0;  P)  =  1: 

er(f;P)  =  Giurt(t)  +  (J  G(t-  a)dU£S)ds 

=  GiUrt(t)  +  C  J  («r  +  K(t  —  S\  P))  ~fi~-ds  (8) 

=  (G+  Qur(t)  +  Giurt(t)  -  t^J  — — — — —  ur(s)ds, 

where  G  =  nrC,  and  with  slightly  more  detail  in  [14].  This  equation  (8)  is  the  general  form  of  the  constitutive  equation 
used  here.  The  value  Go  =  G  +  (  can  be  considered  to  be  a  dynamic  analog  to  the  static  shear  modulus;  this  also 
makes  clear  the  fact  that  Hooke’s  Law  is  incorporated  into  our  model.  We  have  already  discussed  that  Gi  is  the  bulk 
damping  parameter  for  the  Kelvin- Voigt  damping  term.  The  final  integral  represents  a  history  term  which  describes 
the  relaxation  of  the  material  in  response  to  an  applied  stress /strain. 

We  will  ultimately  turn  to  a  discretized  distribution  model  (using  a  discrete  measure  P(r)),  and  connect  it 
to  the  continuum  model  through  a  probability  measure  approximation  as  in  [9].  This  will  allow  us  to  develop  a 
computationally  feasible  inverse  problem,  and  also  give  insight  into  the  underlying  material  mechanics.  But  first  we 
briefly  discuss  a  method  for  approximating  the  loading  process. 
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3.1.5  Approximating  the  loading  process 

Recalling  Figure  2,  the  loading  profile  is  relatively  long  compared  with  the  oscillatory  period;  since  our  concern  is  with 
modeling  the  oscillations,  solving  the  model  from  Ti  is  much  longer  than  necessary.  Also,  early  experimentation  with 
the  model  indicated  that  the  parameters  governing  the  loading  and  resting  process  may  differ  from  those  governing 
the  very  dynamic  post-release  oscillatory  process. 

We  address  these  concerns  by  modeling  the  loading  as  instantaneous  from  at  rest  to  a  displacement  of  A  at 
position  r  =  rmin.  Since  the  material  is  linear,  this  would  then  mean  the  phantom  has  the  profile 


u(r9t) 


A(t max  V'') 
T  max  T  min 


(9) 


up  until  the  time  of  the  weight  release.  Since  this  is  an  approximation,  we  will  neglect  the  times  t  G  (r3,r4),  the 
weight  release  time  period,  since  that  time  frame  is  small  relative  to  the  loading  and  settling  time  from  Ti  to  T3.  We 
also  incorporate  a  time  parameter  T  which  will  represent  our  approximation  of  the  time  when  loading  begins.  In  the 
formulation  here  we  will  use  the  same  relaxation  times  during  the  loading  process  as  during  the  oscillation  period, 
which  means  that  T  has  no  meaning  other  than  as  a  tuning  parameter  that  we  must  estimate.  Thus,  we  assume  the 
given  loading  profiles  for  t  G  (T,  0)  since  T4  =  0  in  our  convention.  This  also  means  that  T  <  0. 

We  incorporate  this  loading  approximation  into  our  model  by  manually  integrating  the  constitutive  relationship 
(8).  For  the  purposes  here,  we  will  call  a  the  full  constitutive  relationship  for  t  >  T  that  is  described  by  (8)  (where 
we  now  use  T  in  the  place  of  Ti),  and  a  the  constitutive  relationship  for  t  >  0.  We  do  this  for  notational  simplicity 
in  the  final  model,  at  the  expense  of  some  minor  notational  confusion  at  the  current  stage.  Note  that  (9)  implies 
that 

ur(r,  t)  = - — - ,  t  G  (T,  0). 

T  max  i min 

Hence,  by  (8)  and  the  above  equation  we  find 

d(t;  P)  =  (G  +  C)ur(t)  +Giurt{t)  -  C  J  d-^L^lUr{s)ds 

^  ^  ^  ^  ,  [°  dK(t  —  s;  P)  ,  ^  f*  dK(t  —  s;  P)  ,  w 

=  (G  +  Q  ur(t)  +  GiUrt(t)  —  C  J  - - ur[s)ds  —  Q  J  - — - ur(s)ds 

=  (G  +  C)ur(t)  +  GlUrt(t)  +  C- - - - {K{t,P)-K{t-T-,P))-C  f  dK{t~S'P)ur(s)ds 

T  max  T min  J  0 

=  a(t;P)-P(f,  T,A,P), 

where  a  and  T  are  respectively  given  by 

cr{t-,P)  =  (G  +  C)  ur(t)  +  Giurt(t)  -  cj  d-^L^Ur{s)ds,  (10) 


and 

p(t;  t,  a,  p )  =  -c — 4 — (*(*; p)  -  -  T; p)- 

T  max  T  min 

By  (10),  we  then  have  the  following: 


®  (J  f  —  (7  f 

•  The  original  stress  boundary  condition  is  (t(rm,;n,  t)  =  0  for  t  >  0.  Using  the  preceding  development,  this 
corresponds  with 

0  =  a(rmin,t ;  P)  -  P{t\  T,  A,  P ) 

which  allows  us  to  write  the  boundary  condition  for  a  model  solved  for  t  >  0  as 


v(rmin,t;  P)  =  P(t;  T,  A,  P). 


We  note  that  the  term 


<7  cr  P(t;T,A,P) 


results  in  a  time-dependent  forcing  term  in  the  shear  PDE. 


r  r  '  r 

We  make  two  comments  before  discussing  the  internal  variable  forms.  First,  if  we  assume,  for  example,  a  single 
relaxation  time  and  that  its  value  is  small,  say  on  the  order  of  1CU1,  then  the  term  K(t  —  T;  P)  =  exp(—  (t—  T)/ri)  ~ 
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exp(— 10(f  —  T))  attains  its  maximum  value  exp(lOT)  when  t  =  0.  Note  that  for,  say,  T  <  — 1,  this  term  is  negligible. 
Relaxation  times  on  this  order  are  what  we  can  later  obtain  in  the  inverse  problem,  which  would  imply  that  in 
our  case  the  material  is  at  rest  after  being  loaded  sufficiently  long  that  it  “forgets”  its  loading  history  by  the  time 
the  weight  is  released.  This  is  good  from  an  experimental  standpoint,  since  the  loading  process  will  never  be  quite 
uniform.  It  is  also  good  to  know  from  a  computational  perspective;  we  can  limit  T  to  being  greater  than  some  value, 
such  as  —20  <  T  <  0,  which  will  keep  the  optimization  algorithm  from  marching  off  unnecessarily  (which  occurred 
in  some  of  our  early  inverse  problem  tests).  Second,  since  we  have  integrated  out  the  loading  history,  we  now  start 
the  model  at  the  time  of  weight  release,  t  =  0.  This  means  that  the  material  is  considered  at  rest  just  prior  to  the 
release;  thus,  in  the  history  integrals  we  will  discuss  in  the  next  section,  all  the  history  now  starts  at  t  =  0  since  the 
history  before  that  point  will  be  incorporated  into  the  initial  loading  profile  and  an  initial  stress  condition. 


3.1.6  Internal  variable  formulation 


In  the  previously  noted  work  on  this  stenosis  problem,  the  double  integrals  that  resulted  from  using  the  continuum  of 
relaxation  times  in  the  stress  equation  were  computationally  intractable  so  another  approach  was  required.  The  idea 
was  to  use  a  discrete  number  of  internal  variables.  As  will  be  noted,  these  gave  rise  to  a  differential  form  which  was 
an  improvement  computationally  since  it  led  to  purely  differential  equations  in  the  model  rather  than  inclusion  of  an 
integro-differential  equation.  With  the  advances  in  desktop  computation  abilities  since  that  time,  the  integral  form 
is  now  reasonable  to  use  in  a  dynamic  model.  However,  internal  variables  are  still  attractive  in  that  they  provide 
a  formulation  that  indicates  some  of  the  internal  material  dynamics.  Physically,  if  we  assume  that  the  molecules 
within  the  biological  tissue  are  on  a  microscopic  scale  then  the  portion  of  the  material  which  is  represented  by  each 
internal  variable  or  internal  strain  ej  is  being  driven  by  the  overall  strain  and  has  a  response  that  varies  depending 
on  the  value  of  the  corresponding  relaxation  time  Tj . 

Previous  work  [8,  42,  50]  assumed  a  discrete  internal  variable  form  as  an  approximation  to  the  Fung  kernel,  using 
the  nonlinear  constitutive  equation  (2).  This  discrete  form  assumed  the  kernel  Q{t)  was  in  an  exponential  form,  with 
the  effects  brought  together  as  a  discrete  sum  in  the  constitutive  equation.  The  results  in  [8,  42,  50]  demonstrate  that 
the  internal  variable  approach  is  valid  and  does  appear  to  work  as  well  as  the  continuum  of  times  in  the  Fung  kernel. 
A  connection  between  the  Fung’s  kernel  and  the  discrete  kernel  is  provided  by  the  work  in  [20].  The  authors  there 

form  the  kernel  Q(t)  =  J  q(t ;  r)dP(r)  where  T  =  [ri,  r2 ]  C  (0,  oo)  is  the  set  of  admissible  relaxation  times,  P(r)  is 

a  probability  measure  on  T ,  and  q(t;  t)  is  a  continuous  function  of  relaxation  times.  If  we  take  q[t ;  r)  =  exp(— t/r), 
this  corresponds  with  the  kernels  previously  discussed.  The  authors  showed  existence  and  uniqueness  results  for  this 
kernel  in  the  nonlinear  constitutive  equation  (2).  Then,  a  result  from  [9]  allows  one  to  approximate  any  measure  P(r) 
with  a  discrete  measure.  This  discrete  measure  approximation  leads  us  back  to  the  case  with  a  sum  of  exponentials, 
but  from  the  probabilistic  framework  we  know  conclusively  that  we  are  approximating  the  continuous  spectrum  of 
Fung’s  kernel  and  that  this  approximation  has  been  viable  when  implemented. 

With  this  background  on  previous  work  using  internal  variables  in  hand,  we  move  forward  by  modifying  our 
current  model.  We  manipulate  the  form  of  Equation  (10)  as  follows: 


cr(t;P) 


(G  +  C)  Ur(t)  +  G\Urt(t)  -  C  [  dK{\  S'P)Ur{s)ds 

Jo  9s 

(G  +  C)  Ur{t)  +  Gmnit)  -  C  J  exv(-(t  -  s)/t)gLP(t)^  Ur(s)ds 

(G  +  C)  ur(t)  +Giurt(t)  -  C  J  £i (t;r)dP(r), 


(11) 


where  we  let  ei(t;r) 
differential  form 


r  d 

/  —  (exp  (— (t  —  s)/t))  ur(s)ds.  Rather  than  the  integral  form  for  e±,  we  can  use  the 

Jo  9s 


Tdtei^’T)  +  ei^;r) 


ur(t),  €i(0;t)=0 


(12) 


which  is  then  solved  simultaneously  with  the  rest  of  the  model  dynamics.  This  is  then  an  internal  variable  or  internal 
strain ,  driven  by  the  overall  strain  ur{t)1  which  is  the  continuous  form  of  the  internal  variable  formulation. 

We  now  may  finally  make  the  discrete  assumption 


Np 

\  " 


P(T)  =  }^PiA 


i= i 
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where  AT.  is  the  Heaviside  function  with  step  at  Tj ,  and  pj  are  the  proportions  of  the  material  subject  to  relaxation 

Np 

time  Tj  so  that  YY^i  =  1.  By  substituting  this  discrete  P  into  the  form  for  a  as  developed  in  (11),  we  obtain  the 

i= i 

discrete,  internal  variable  form  of  the  constitutive  relationship 


/  Np  \  Np 

a(t)  =  G  +  Ur{t)  +  Gmrtit)  -  ^(jp (t), 

V  j= i  / 


(13a) 


with  internal  variables  obeying  (for  j  =  1,2, ,  Np ) 

Tj^e3(t)  +  e:l{t)  =  Ur(t),  eJ(0)  =  0,  (13b) 

Np 

and  where  we  have  defined  Cj  =  C Pj  so  that  £  =  YY  Q,  and  P  =  t\ (•;  Tj). 

1=1 


3.2  Final  shear  wave  model 

We  now  put  together  the  shear  wave  equation  (1),  using  the  constitutive  equation  (13)  but  with  the  loading  history 
approximation  incorporated  as  discussed  in  Section  3.1.5.  Recall  also  that  the  discrete  assumption  for  P  and  the 
form  of  K  gives  us 

C K{t ;  P)  =  C  T,f=i  Pj  exP i-t/Tj)  =  Cl  exp(— t/rj). 

These  equations  are  just  manipulated  versions  of  the  general  equations  of  Theorem  3.1,  so  we  still  know  a  unique 
weak  solution  exists  on  any  finite  time  interval. 

The  shear  wave  equation,  solved  for  t  >  0  which  is  the  release  time,  are  then 


d2  d  cr(r,t)  1  A 

PT> 72u(r’t)  ~  Ta(r’t) - = - 

otz  dr  r  r  rmax  -  rmm 


Np 


Np 


^iXmin  i  t)  — 


-A 


Tmax  ' min 


Np  Np 

Cl  exp  i-t/Tj)  -  Y2  (j  exp(— (t  -  r)/Tj) 
l=i  l=i 


Ci  exp(— t/ri)  -  Y2  Ci  exp  -  T)/ri) 

=1 

,  u(rmax,t)  =  0 


i  (-.n,  A(rmax  r )  .  . 

u{r,  0)  = - ,  wt(r:0)=°i 


T  max  T  min 


where 

Np  \  Np 

G  +  Cl  J  ur(t)  +  G\urt{t)  —  YY  (t) 
1=1  /  1=1 

with  the  internal  variables  subject  to  (for  j  =  1,2,...,  Np) 

Tj^3(t)  +  £:’(t)=Ur,  eJ(0)  =  0. 


(14a) 


(14b) 


(14c) 


Np 

We  note  that  Go  =  G  +  YY  Cl  is  the  dynamic  analog  of  the  shear  modulus. 

1=1 

3.3  Numerical  method 

The  shear  wave  model  (14)  is  numerically  solved  by  using  the  high  order  space-time  finite  element  method.  Specifi¬ 
cally,  in  time,  we  use  a  discontinuous  Galerkin  method  composed  of  normalized  Legendre  polynomials  (of  order  4).  In 
space,  we  use  a  continuous  spectral  finite  element  method  composed  of  Lagrange  basis  functions  on  Gauss-Lobatto 
nodes  (also  of  order  4).  Under  this  scheme,  the  system  matrices  are  diagonalizable  (this  could  be  lost  if  we  did 
not  use  normalized  Legendre  polynomials  in  time),  and  hence  the  time-coupled  computations  within  a  time  step 
can  be  decoupled.  This  makes  the  reasonably  high  order  finite  element  time  discretizations  feasible.  An  argument 
for  adopting  such  higher  order  schemes  instead  of  the  lower  order  ones  is  that  higher  order  spatial  discretizations 
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have  been  found  highly  desirable  for  wave  equations  in  terms  of  the  control  of  dispersion  errors  (e.g.,  see  [1,  23]). 
Moreover,  higher  order  schemes  are  capable  of  providing  higher  fidelity  solutions  than  lower  order  schemes  for  the 
same  amount  of  computational  work.  This  is  especially  important  to  make  the  inverse  problem  practical  (as  each 
inverse  problem  may  require  solving  the  forward  problem  (i.e.,  shear  wave  model)  hundreds  of  times).  Further  details 
on  this  numerical  method  are  in  a  forthcoming  BICOM  report  [38],  wherein  an  extensive  set  of  numerical  results 
are  given  to  demonstrate  the  favorable  effect  on  the  numerical  error  and  computational  work  of  the  higher  order 
temporal  and  spatio-temporal  approximations. 


4  Inverse  Problem 

With  models  in  hand,  we  now  turn  to  matching  the  model  output  to  data.  We  will  use  two  common  methods  in 
order  to  estimate  model  parameters.  One  is  ordinary  least  squares  (OLS)  and  the  other  is  generalized  least  squares 
(GLS).  These  will  be  defined  later  in  Section  4.1. 

As  discussed  in  Section  2,  an  experiment  has  been  designed  to  gather  one-dimensional  shear  data.  Measurements 
in  our  experiment  are  taken  at  r  =  rmin,  and  will  be  denoted  by  Uj  at  measurement  time  point  tj.  Corresponding 
shear  model  solutions  at  the  same  spatial  location  will  be  denoted  u{tj ;  10e),  where  the  measurement  location  has  been 
suppressed  for  notation  convenience  and  where  9  represents  a  vector  of  the  base-10  logarithm  of  each  parameter  (the 
same  idea  used  previously  [13]  to  reduce  parameter  scaling  issues).  The  data  set  has  been  trimmed  to  the  dynamic 
oscillations  after  the  release,  and  thus  the  time  frame  is  200  ms.  The  data  were  sampled  at  a  rate  of  2048  Hz; 
however,  this  high  rate  proved  to  make  the  inverse  problem  difficult  and  computationally  intractable  because  that 
many  data  points  resulted  in  the  inverse  problem  being  over  determined.  Instead,  we  will  use  every  other  data  point 
from  the  larger  data  set  for  a  sampling  rate  of  1024  Hz  which  will  we  later  refer  to  as  the  “every  data  point”  set. 
We  take  n  to  be  the  total  number  of  data  points  for  a  particular  data  set,  and  thus  can  describe  the  measurement 
time  points  for  the  full  “every  data  point”  set  as  tj  =  j /1024  where  j  =  0, 1, . . . ,  n  —  1.  There  will  also  be  a  reduced 
data  set  where  we  take  every  other  data  point  starting  with  to  =  0;  which  corresponds  with  a  data  sampling  rate  of 
512Hz. 

Since  some  of  the  data  points  were  near  zero  in  absolute  value,  we  found  that  those  points  resulted  in  scaling 
problems  when  using  the  GLS  model  to  estimate  model  parameters  (since  the  corresponding  cost  functional  divides  by 
the  model  value  as  we  will  see  later  when  this  method  is  defined).  To  account  for  this,  we  removed  from  consideration 
any  data  points  Uj  where  \uj\  <  5  x  10~6.  This  value  was  chosen  by  examining  the  data,  noting  that  the  data  is  on 
the  order  of  10-5  and  that  the  “jitter”  one  can  see  in  Figure  2  has  a  magnitude  of  roughly  5  x  ICG6  during  the  times 
before  loading  up  to  Ti,  then  during  the  settling  period  from  T2  to  T3,  and  again  in  the  settling  period  after  the 
oscillations  have  died  out.  Thus,  our  threshold  level  is  below  the  level  of  noise  in  the  data.  This  level  also  eliminated 
only  a  few  data  points,  while  providing  significantly  improved  GLS  robustness.  The  number  of  data  points  n  is  then 
reduced  according  to  how  many  thresholded  data  points  were  removed. 

Before  going  into  the  setup  and  results  for  the  inverse  problem,  we  note  that  the  forward  (i.e.,  direct)  problem 
where  we  solve  for  displacement  (using  the  method  discussed  in  Section  3.3)  is  as  follows: 

•  Forward  problem:  Given  G,  G i,  Tj  and  ( j  for  j  =  1,2, ,  Np.  T,  A,  rmin,  rmax ,  and  p ,  solve  model  (14)  for 
displacement  u{r,t)  at  each  position  x  £  [rmax,rmjn]  for  t  £  [0,T]. 

The  inverse  problem  we  will  develop  here  is  as  follows: 

•  Inverse  problem:  Given  shear  displacement  data  at  r  =  rmin  and  a  corresponding  forward  problem  solver 
for  displacement,  along  with  specified  values  for  p,  rmin,  and  rmax,  find  values  for  the  constants  G,  G i,  Tj  and 
Q  (for  j  =  1,2,...,  Np),  A,  and  T  which  provide  the  best  fit  to  the  data  (in  a  manner  which  will  be  defined 
shortly) . 

We  assume  that  the  parameters  lie  in  some  admissible  set  Q  C  KK,  where  Q  is  assumed  to  be  compact  and  k  is 
the  number  of  parameters  requiring  estimation.  Throughout  the  remainder  of  this  work,  we  will  denote  the  log-scaled 
parameter  vector  (for  Np  =  1  and  k  =  6)  as 

9  =  (log10(G),log10(Gi),log10(Ci),log10(n),log10(-A),log10(-T)).  (15) 

Thus,  as  long  as  we  define  our  cost  function  to  be  a  continuous  function  of  the  parameters,  we  know  the  inverse 
problem  has  a  solution  (minimizing  a  continuous  function  on  a  compact  parameter  space).  One  could  broaden  this 
parameter  estimation  framework  to  the  distributional  case  if  desired,  taking  an  admissible  parameter  space  as  a 
compact  subset  of  Euclidean  space  (including  all  parameters  excuding  relaxation  times)  along  with  the  space  of 
probability  measures,  and  use  the  Prohorov  metric  framework  (see,  e.g.,  [7,  15,  Sec.  4])  and  the  approximation 
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results  of  [9].  This  again  leads  to  minimizing  a  continuous  function  of  the  parameters  over  a  compact  space.  Either 
way,  the  inverse  problems  we  will  shortly  define  will  have  solutions. 

4.1  Statistical  Models  and  Parameter  Estimators 

In  order  to  carefully  define  the  way  in  which  we  will  measure  the  closeness  of  the  data  to  model  values,  we  must  first 
discuss  underlying  statistical  models  for  the  error  present  in  the  data.  A  proper  error  model  is  also  key  to  correctly 
determining  parameter  confidence  intervals.  Much  of  the  discussion  here  is  similar  to  that  in  [13],  with  background 
on  ordinary  least  squares  (OLS)  and  generalized/ weighted  least  squares  (GLS  or  WLS)  given  in  [22],  for  example. 

We  will  assume  the  errors  £j  are  independent,  identically  distributed  with  mean  zero  {E\£:)]  =  0)  and  constant 
variance  var (£j)  =  a 2;  this  process  has  realizations  £j.  Note  that  we  do  not  assume  we  know  the  underlying 
distributions  from  which  the  errors  come;  we  only  know  the  first  two  central  moments  as  specified.  We  use  this  error 
process  in  proposing  two  error  models  and  corresponding  parameter  estimators. 

•  Absolute  error:  Here  we  have  the  error  process  Uj  =  u(tj\  lO6*0)  +  £:r  with  realizations 

Uj  =  u(tj\  109°)  +  £j,  (16) 

where  90  is  some  hypothesized  “true”  parameter  value  (see  [22]).  We  use  the  ordinary  least  squares  cost  function 

n— 1 

Jols(9)  =  J2luj~u(tE  109)]2- 

3=  0 

The  corresponding  inverse  problem  for  the  logged  parameters  is  then 

n—  1 

Oois  =  argmin  Jou(0)  =  argmin  Y,iu3  -  u(tj;  10s)]2.  (17) 

8eQ  eeQ  w 

3=0 

This  function  minimizes  the  distance  between  the  data  and  model  where  all  observations  are  considered  to  have 
equal  importance  (weight).  Since  u(tj\  10s)  is  a  continuous  function  of  9 ,  J0\s  is  also  a  continuous  function 
of  9 ,  which  means  we  are  minimizing  a  continuous  function  of  9  over  a  compact  set  Q ,  and  thus  this  inverse 
problem  has  a  solution. 

•  Relative  error:  Here  we  have  the  error  process  Uj  =  u(tj ;  10s°)  +  u(tj ;  10 S°)£j  with  realizations 

Uj  =  u(tj-  10s°)  +  u(tj ;  10  6o)sj-  (18) 

For  this  case,  we  construct  the  generalized  (weighted)  least  squares  cost  function  (as  per,  e.g.,  [22]) 

n—  1 

Jgis(o)  =  Y w2jluj  - 1Qe)]2 

3=0 

where  we  define  the  weights  Wj  =  u(tj\  10s)-1.  In  this  case,  since  we  are  examining  a  relative  error  model 
(18),  these  weights  take  into  account  the  unequal  quality  of  observations;  dividing  by  the  function  value  has  a 
“normalizing”  effect  on  the  errors,  accounting  for  the  scale  differences  which  may  be  present  in  the  errors  at 
larger  versus  smaller  model  values. 

We  now  wish  to  find  9  such  that  £fgis{9)  is  minimized.  We  can  either  solve  this  directly,  or  by  using  an  iterative 
procedure  in  order  to  estimate  9gis  (since  the  weights  must  also  be  estimated).  We  will  use  an  iterative  method, 
described  as  follows  (see  [22]  and  references  therein  for  convergence  details): 

1.  Define  9°  =  9ais,  and  set  k  =  0. 

2.  Form  the  weights  uij  =  u(tj ;  10s  )-1,  using  weight  thresholding  (described  below). 

3.  Re-estimate  9gis  by  solving 

n—1 

9k+1  =  argmin  Y  ~  u(Ab  10<?)]2 

W  3=0 

to  obtain  the  k  +  1  estimate  9k+1  for  9gis. 
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4.  Set  k  =  k  +  1  and  return  to  Step  2.  Terminate  when  successive  estimates  for  9gis  are  sufficiently  close, 
or  when  one  has  iterated  20  times.  For  our  problem,  our  “sufficiently  close”  criterion  was  found  by 
determining  if  ||0fe+1  —  0fc||oo  <  10~3,  where  ||0||oo  is  the  maximum  component  of  the  given  vector  9.  The 
parameter  values  being  estimated  are  all  log-scaled,  and  are  thus  on  the  order  of  [10— 1 , 101] .  This  puts 
the  stopping  criterion  at  two  orders  of  magnitude  less  than  the  smallest  log-scaled  parameter  value,  which 
is  sufficient  in  our  problem. 

Even  though  we  have  removed  all  data  points  with  absolute  value  under  5  x  10~6,  we  still  account  for  the  (now 
unlikely)  possibility  that  some  model  values  may  still  end  up  small  in  absolute  value.  Thus,  we  incorporate 
thresholding  on  the  weights  to  keep  from  dividing  by  zero.  We  take  a  threshold  value  of  1  x  10-1°,  as  this  is 
almost  certainly  below  the  threshold  of  significance  in  terms  of  the  model  displacements.  Then,  for  all  indices 
j  €  {k  |  | u(tk',  10e)|  <  1  x  10~10},  we  set  wj  =  1  x  1010.  This  is  done  each  time  the  weights  are  re-estimated 
in  Step  2  of  the  iterative  process. 

With  weight  thresholding,  we  are  assured  that  the  iterative  process  is  possible  numerically.  Thus,  similar  to 
the  ordinary  least  squares  case,  at  each  step  k  in  the  iterative  GLS  estimation  process  we  are  minimizing  a 
continuous  function  of  9  over  a  compact  parameter  space  Q ,  and  thus  the  inverse  problem  in  each  iteration 
will  have  a  solution.  Also,  as  long  as  the  iterative  process  is  carried  out  sufficiently  many  times,  under  certain 
conditions  the  weights  will  converge  uij  — >•  u(tg\  lO6*®'8)-1  (see,  e.g.,  [22]). 

4.1.1  Optimization  considerations 

As  in  [13],  we  used  the  Matlab  routine  lsqnonlin  for  our  optimization  routine  to  solve  for  90is  and  9gis.  We  used  the 
trust-region-reflective  (TRR)  algorithm  that  is  built  in;  as  our  previous  effort  in  [13]  demonstrated,  the  Levenburg- 
Marquardt  option  was  slower  than  TRR  and  did  not  give  us  better  results.  Since  we  are  using  at  least  one  relaxation 
time,  we  do  not  consider  fmincon  which  we  have  shown  to  be  ineffective  in  estimating  relaxation  times. 

In  order  to  start  the  optimization  routines  for  computing  90is ,  we  must  provide  initial  parameter  values  (for  9gis 
we  use  the  estimated  value  for  90is  as  our  initial  guess).  From  a  perusal  of  the  viscoelastic  materials  literature,  our 
experience  with  the  previous  conceptual  work,  and  from  some  manual  examination  on  the  current  data  sets,  we  chose 
the  initial  values  as  follows: 

G  =  4.5  x  103  Pa,  Gi  =  5  Pa  •  s,  Ci  =  2.8  x  104  Pa,  n  =  0.06  s,  A  =  -1.7  x  1(T4  m,  T  =  -0.01  s. 

As  log-scaled  values  (c.f.  (15)),  this  gives  us 

6°gls  =  (3.6532,0.6990,4.4472,-1.2218,-3.7696, -2)t. 


4.1.2  Residuals 

We  will  also  include  residual  plots  to  assist  in  analysis  of  the  model  fit  to  data,  and  to  indicate  which  error  model 
best  describes  the  error  in  the  data.  Residuals  give  a  sense  for  the  model  fit  to  data,  but  more  importantly  the 
residuals  can  give  an  indication  [22]  regarding  the  appropriateness  of  our  error  model.  If  the  absolute  residuals  seem 
to  be  randomly  dispersed  around  the  horizontal  axis  and  form  a  horizontal  band  around  that  axis,  then  the  absolute 
error  model  may  be  correct.  On  the  other  hand,  if  the  (modified)  relative  residuals  seem  to  be  randomly  dispersed, 
then  the  relative  error  model  may  be  correct.  We  define  the  following: 

•  Absolute  residuals  are  computed  as  r;-  =  Uj  —  u(tj\  106*),  where  9  is  the  particular  parameter  estimate  being 
considered. 

•  Relative  residuals  are  computed  as  rg  =  Wj(uj  —  u(tj ;  10s))  where  Wj  =  u(tj ;  106*)-1  and  the  Wj  are  thresh- 
olded  in  the  same  manner  as  discussed  earlier. 

4.2  Inverse  problem  results,  Np  —  1 

We  now  demonstrate  the  ability  of  our  model  to  match  data.  For  this  purpose,  we  will  take  a  single  relaxation  time 
(Np  =  1)  as  that  is  enough  to  show  model  fidelity  to  data.  We  run  both  the  absolute  (OLS)  and  relative  (GLS)  error 
models  on  a  sample  data  set  using  a  264  g  loading  weight.  We  will  report  parameter  estimates,  standard  errors  and 
confidence  intervals,  plots  of  model  fits  to  data,  plots  of  residuals  versus  time,  and  plots  of  residuals  versus  model 
values.  We  use  these  elements  in  order  to  recommend  error  models. 
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Standard  errors  (and  corresponding  confidence  intervals)  are  computed  using  asymptotic  error  theory.  For  the 
absolute  error  model,  the  process  is  the  same  as  that  which  we  used  in  [13],  and  is  also  described  in  [22,  Ch.  3]; 
for  the  relative  error  model,  the  corresponding  asymptotic  error  methodology  is  discussed  in  [22,  Ch.  3].  Since  the 
theory  is  common  enough,  we  do  not  reiterate  it  here  and  refer  interested  readers  to  the  aforementioned  references. 

We  have  also  examined  parameter  estimation  using  data  sampled  at  different  rates.  This  allowed  for  a  study  of 
whether  the  parameter  estimates  and  associated  confidence  intervals  remain  consistent  as  the  number  of  data  points 
is  reduced.  Using  fewer  data  points  is  also  a  way  of  decreasing  computational  times  for  the  inverse  problems.  We 
ran  the  inverse  problem  on  each  data  set  and  using  each  error  model  using  all  the  data  points  (1024  Hz)  and  using 
every  other  data  point  (512  Hz),  as  discussed  at  the  beginning  of  Section  4. 

Numerical  results  show  that  the  results  for  data  sampled  at  1024  Hz  is  consistent  with  those  found  for  the  data 
sampled  at  512  Hz.  Hence,  we  report  only  the  results  for  data  sampled  at  512  Hz  (all  results  for  data  sampling 
at  1024  Hz  are  available  in  [14]).  The  parameter  estimates  and  confidence  intervals  (see  [10,  11,  13,  16,  22]  for 
information  on  computing  confidence  intervals),  obtained  from  the  inverse  problems  using  OLS  and  GLS,  are  shown 
in  Tables  1-2.  We  observe  from  these  two  tables  that  the  standard  errors  for  the  GLS  case  is  larger  than  in  the 


Table  1:  Shear  optimization  results  and  confidence  analysis  for  OLS  on  a  264  g  data  set  using  every  other  data  point. 


Param. 

Estimate 

SE 

CI95 

iogioCG) 

3.5431 

0.2498 

(3.0469,  4.0393) 

logio(Gi) 

0.3753 

0.3227 

(-0.2657,  1.0164) 

l°gio(Ti) 

-1.4450 

0.2474 

(-1.9364,  -0.9535) 

!ogio(Ci) 

4.4761 

0.0294 

(4.4178,  4.5345) 

l°gio  {-A) 

-3.7501 

0.0071 

(-3.7643,  -3.7360) 

!ogi0(-T) 

-2.1813 

0.2779 

(-2.7334,  -1.6293) 

Shear  modulus  dynamic  analog  Go  =  33.423  kPa 


Table  2:  Shear  optimization  results  and  confidence  analysis  for  GLS  on  a  264  g  data  set  using  every  other  data  point. 


Param. 

Estimate 

SE 

CI95 

!ogio(G) 

3.8649 

1.0435 

(1.7922,  5.9377) 

logio(Gi) 

0.5449 

0.3561 

(-0.1625,  1.2523) 

l°gio(ri) 

-1.0217 

1.0543 

(-3.1161,  1.0726) 

l°gio(Ci) 

4.4171 

0.2918 

(3.8375,  4.9967) 

l°gio(— A) 

-3.8026 

0.0119 

(-3.8263,  -3.7789) 

!ogi0(-T) 

-1.6729 

1.3806 

(-4.4153,  1.0695) 

Shear  modulus  dynamic  analog  Go  =  33.454  kPa 


OLS  case.  In  addition,  we  see  that  the  standard  error  for  G\  is  on  the  same  order  of  magnitude  as  the  parameter 
estimate  itself  for  both  the  OLS  and  GLS  results,  and  the  standard  errors  for  t\  and  T  in  the  GLS  case  are  also  on 
the  same  order  of  magnitude  as  their  corresponding  estimates.  This  is  consistent  with  the  sensitivity  results,  where 
the  model  output  is  less  sensitive  to  G i,  t\  and  T  than  to  G,  £i,  and  A  (more  details  are  available  in  [14]). 

Model  fits  as  well  as  residual  plots  are  shown  in  Figure  3.  In  all  cases,  the  model  fits  to  data  are  good.  In  the 
bottom  row  of  Figure  3,  the  residuals  versus  model  plots  are  not  noticeably  different  between  the  OLS  and  GLS 
cases.  The  initial  indication  is  that  we  have  more  confidence  in  the  OLS  results.  However,  the  residual  versus  time 
plots  (middle  row  of  Figure  3)  raise  cause  for  concern.  In  the  OLS  residual  versus  time  plots,  there  is  a  noticeable 
“fan”  structure  for  early  times.  However,  for  the  GLS  error  model,  the  residual  versus  time  plots  do  not  show  a  fan 
structure  and  are  fairly  randomly  distributed.  Since  this  indicates  that  the  OLS  error  model  may  not  be  correct, 
we  are  inclined  to  recommend  the  GLS  error  model  so  that  we  do  not  mistakenly  overstate  our  confidence  in  the 
parameter  estimates,  which  we  could  do  if  we  used  the  parameter  estimates  from  the  possibly-wrong  OLS  case. 
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Abs.  residuals  (units:  m)  Abs.  residuals  (units:  m)  Displacement  (units:  m) 


Shear  optimized  fit,  OLS,  N  =1  (264g,  trial  3) 
x  10  p 


Shear  residuals  vs  time,  OLS,  N  =1  (264g,  trial  3) 
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Shear  optimized  fit,  GLS,  N  =1  (264g,  trial  3) 
x  10  p 


Shear  residuals  vs  time,  GLS,  Np=1  (264g,  trial  3) 
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Shear  residuals  vs  model,  OLS,  N  =1  (264g,  trial  3) 
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Figure  3:  Shear  data  fit  using  data  sampled  at  512  Hz,  Np  =  1,  weight  264  g.  (left)  OLS  results. 


(right)  GLS  results. 
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5  Discussion  and  Future  Work 


We  have  developed  an  updated  one-dimensional  viscoelastic  model  for  tissue  and  have  used  experimental  data  from  a 
simple  homogeneous  gel  phantom  to  test  the  ability  of  our  model  to  describe  wave  propagation  in  the  medium.  The 
data  were  generated  from  a  drop  experiment  designed  to  produce  oscillations  in  the  gel  of  a  magnitude  comparable 
to  that  produced  by  blood  flow  in  a  stenosed  coronary  artery  impacting  the  vessel  wall,  a  disturbance  which  results 
in  shear  wave  propagating  away  from  the  vessel  walls  downstream  of  the  blockage.  In  our  inverse  problem  results  as 
discussed  in  Section  4.2,  we  have  shown  an  ability  to  consistently  model  the  wave  propagation  using  different  error 
models  and  at  different  data  sampling  frequencies,  obtaining  good  fits  to  data  in  all  of  our  inverse  problems.  In 
addition  to  a  good  fit,  though,  we  also  examined  statistical  properties  of  the  parameter  estimators  as  well  as  residual 
plots  to  gain  more  insight  into  the  proper  error  model  for  the  shear  data  set.  This  is  necessary,  since  a  correct  error 
model  is  essential  in  order  to  apply  the  asymptotic  error  theory  properly  and  thus  obtain  correct  confidence  intervals. 
We  recommend  taking  the  more  conservative  route  and  using  the  GLS  parameter  estimates;  even  though  the  GLS 
estimates  had  larger  standard  errors,  there  were  indications  from  the  residual  versus  time  plots  for  OLS  that  the  OLS 
model  is  not  correct.  Overall,  we  have  successfully  demonstrated  the  ability  of  mathematical  model  to  accurately 
describe  the  data  from  laboratory  experiments  using  a  homogenous  tissue-mimicking  material  gel  phantom.  A  linear 
viscoelastic  constitutive  relationship,  i.e.,  (14b),  was  adequate.  This  is  a  significant  achievement,  as  all  the  work 
previously  discussed  was  limited  to  inverse  problems  on  simulated  data  or  data  that  was  not  from  the  impulse-type 
experiments. 

We  do  note  that  both  the  model  and  experiments  that  we  studied  here,  though  quite  useful  as  a  starting  point, 
are  limited  in  replicating  reality.  Two  different  paths  are  the  likely  next  steps.  First,  we  are  currently  examining 
a  two-dimensional  model  and  corresponding  experimental  configurations.  Experiments  are  currently  in  progress  to 
produce  a  two-dimensional  wave  from  different  points  in  the  medium  and  with  different  detection  points  along  the 
outer  wall  of  the  phantom.  It  is  conceivable  that  the  one-dimensional  parameters  could  be  used  as  a  rough  first 
approximation  in  a  corresponding  two-dimensional  code,  which  would  allow  us  to  focus  on  trying  to  determine  the 
location  of  the  wave  generation  in  the  medium.  Also,  these  parameter  values  could  be  used  in  a  model  of  wave 
propagation  in  another  conceptual  device  designed  to  mimic  a  constricted  artery  and  the  waves  that  result  from 
passing  fluid  through  a  constricted  pipe  in  the  center  of  the  medium.  Therefore  in  the  slightly  longer  term,  we 
will  also  likely  need  to  conduct  an  inverse  problem  using  a  two-dimensional  model  and  corresponding  data.  These 
one-dimensional  results  will  provide  a  starting  point  for  parameters  in  that  inverse  problem,  hopefully  decreasing 
runtime  and  the  time  it  takes  to  find  viable  parameters.  The  same  issues  discussed  here  (sensitivity  to  parameters, 
data  frequency,  number  of  relaxation  times)  will  again  be  of  concern  for  the  two-dimensional  problem.  Future  efforts 
will  also  involve  scaling  up  all  these  experiments  to  larger  phantoms  and  then  to  some  sort  of  actual  tissue  sample 
experiments. 

Second,  the  work  to  this  point  has  still  been  focused  on  the  simplest  problem,  that  of  a  homogeneous  medium.  A 
clear  next  step  would  be  adding  inhomogeneous  features  to  the  gel  phantom  to  mimic  different  types  of  tissue  in  the 
chest  cavity.  This  would  likely  be  formulated  initially  in  the  one-dimensional  case,  so  that  we  could  study  the  model 
responses  in  a  simpler  framework  before  moving  to  a  more  complicated  two-dimensional  case.  We  can  increase  the 
number  of  internal  variables  to  accommodate  different  material  relaxation  properties.  The  model  also  allows  changes 
to  the  material  parameters;  for  example,  we  could  use  piecewise-defined  constants  to  represent  different  types  of 
tissue.  These  two  options  represent  reasonable  directions  for  the  next  stages  of  investigation. 
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